load("./data/elk_processed.rda")
library(sf)
library(dplyr)
elk_res <- elk_processed |>
filter(id %in% c("YL58","YL64","YL80", "YL91", "YL94"))
elk_res_sf <- elk_res |>
st_as_sf(coords = c("lon","lat"), crs=4326) |>
st_transform(32611)
trails_and_roads <- st_read("./Hebblewhitematerials/data/humanacess.shp") |>
st_transform(32611)
## Reading layer `humanacess' from data source
## `C:\Users\nicol\Documents\BrazilMove2024\Hebblewhitematerials\data\humanacess.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 24134 features and 39 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 513240.9 ymin: 5665028 xmax: 650430.4 ymax: 5760389
## Projected CRS: NAD83 / UTM zone 11N
str(trails_and_roads)
## Classes 'sf' and 'data.frame': 24134 obs. of 40 variables:
## $ IGDS_LAYER: chr NA NA NA NA ...
## $ IGDS_LEVEL: int 0 0 0 0 0 0 0 0 0 0 ...
## $ IGDS_COLOR: int 0 0 0 0 0 0 0 0 0 0 ...
## $ NAME : chr NA NA "Spreading/ Porcupine conn" "Porcupine Cr. Tr." ...
## $ CODE : chr "DC31700100" "DC31700100" "DC31700100" "DC31700100" ...
## $ TYPE : chr "TRAIL" "Trail" "Trail" "Trail" ...
## $ CLASS : chr NA NA NA NA ...
## $ SOURCE : chr "BAG" "OldHUM" "OldHUM" "OldHUM" ...
## $ ACCURACY : num 0 0 0 0 0 0 0 0 0 0 ...
## $ STATUS : chr NA NA NA NA ...
## $ NOTES : chr NA NA NA NA ...
## $ MOTOR : int 0 0 0 0 0 0 0 0 0 0 ...
## $ OHV : int 0 0 0 0 0 0 0 0 0 0 ...
## $ FOOT : int 1 1 1 1 1 1 1 1 1 1 ...
## $ CYCLE : int 0 0 0 0 0 0 0 0 1 1 ...
## $ HORSE : int 0 0 0 0 0 0 1 1 1 1 ...
## $ SKI : int 0 0 0 0 0 0 0 0 1 1 ...
## $ SNOWM : int 0 0 0 0 0 0 0 0 0 0 ...
## $ COMMERCIAL: int 0 0 0 0 0 0 0 0 0 0 ...
## $ OTHER : chr NA NA NA NA ...
## $ SUM_CLASS : chr "MEDIUM" "MEDIUM" "MEDIUM" "MEDIUM" ...
## $ SUM_USE : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ESTIMATE : chr "BAG" "BAG" "BAG" "BAG" ...
## $ CONFIDENCE: num 0 0 0 0 0 0 0 0 0 0 ...
## $ SUM_YEAR : chr "2001" "2002" "2002" "2002" ...
## $ WIN_CLASS : chr "0" NA NA NA ...
## $ WIN_USE : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ESTIMATE0 : chr "0" "0" "0" "0" ...
## $ CONFIDEN0 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ WIN_YEAR : chr "0" "0" "0" "0" ...
## $ MAP_SHEET : chr NA NA NA NA ...
## $ YEAR_MAP : chr NA "2001" "2001" "2001" ...
## $ YEAR_UPD : chr "2002" "2002" "2002" "2002" ...
## $ DISPOSITIO: chr NA NA NA NA ...
## $ SOURCETHM : chr NA NA NA NA ...
## $ MOTO_YEAR : chr NA NA NA NA ...
## $ YEAR : chr NA NA NA NA ...
## $ COMMENT : chr NA NA NA NA ...
## $ LENGTH : num 1013 1469 4149 4079 3988 ...
## $ geometry :sfc_MULTILINESTRING of length 24134; first list element: List of 1
## ..$ : num [1:6, 1:2] 543862 544276 544400 544438 544469 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTILINESTRING" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:39] "IGDS_LAYER" "IGDS_LEVEL" "IGDS_COLOR" "NAME" ...
st_crs(trails_and_roads)==st_crs(elk_res_sf)
## [1] TRUE
plot(st_geometry(trails_and_roads))
plot(elk_res_sf$geometry, add=TRUE, col="green")
library(raster)
landcover <- raster("./Hebblewhitematerials/data/landcover.tif")
landcover
## class : RasterLayer
## dimensions : 3811, 2652, 10106772 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 540827.8, 620387.8, 5642002, 5756332 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
## source : landcover.tif
## names : TRUE.
## values : 0, 16 (min, max)
plot(landcover)
landcover <- projectRaster(landcover, crs = crs(elk_res_sf))
0 - NA 1 - Open Conifer Forest 2 - Moderate Conifer Forest 3 - Closed Conifer Forest 4 - Deciduous Forest 5 - Mixed Forest 6 - Regeneration 7 - Herbaceous 8 - Shrub 9 - Water 10 - Rock-Ice 11 - Cloud 12 - Burn-Forest 13 - Burn-Grassland 14 - Burn-Shrub 15 - Alpine Herb 16 - Alpine Shrub
sort(unique(round(landcover@data@values)))
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
landcover_rc <- reclassify(landcover, c(-Inf, 0, 0,
0, 1, 1,
1, 2, 1,
2, 3, 1,
3, 4, 1,
4, 5, 1,
5, 6, 2,
6, 7, 3,
7, 8, 3,
8, 9, 4,
9, 10, 5,
10, 11, 6,
11, 12, 7,
12, 13, 7,
13, 14, 7,
14, 15, 8,
15, Inf, 8))
0 - NA 1 - Forest (Conifer, Deciduous, Mixed) 2 - Regeneration 3 - Herbaceous/Shrub 4 - Water 5 - Rock-Ice 6 - Cloud 7 - Burned (Forest, Grassland, Shrub) 8 - Alpine (Herb, Shrub)
sort(unique(round(landcover_rc@data@values)))
## [1] 0 1 2 3 4 5 6 7 8
plot(landcover_rc)
plot(st_geometry(elk_res_sf), add=TRUE)
library(sp)
library(adehabitatHR)
elk_mcp <- mcp(as_Spatial(elk_res_sf), 100)
determineAvailableLocations <- function(used_data, mcp){
random_n <- nrow(used_data)+1000
sample <- spsample(mcp, random_n, "random") |> st_as_sf()
sample$id <- used_data$id[1]
sample$Used <- FALSE
used_data$Used <- TRUE
sample <- sample |> dplyr::select(id, Used, geometry)
used_data <- used_data |> dplyr::select(id, Used, geometry)
combine_data <- rbind(used_data, sample)
return(combine_data)
}
elk_id <- split(elk_res_sf, elk_res_sf$id)
elk_id_combine <- lapply(elk_id,
determineAvailableLocations,
elk_mcp)
elk_combine <- do.call("rbind", elk_id_combine)
elk_mcp_sf <- elk_mcp |>
st_as_sf()
trails_and_roads2 <- st_crop(trails_and_roads, elk_mcp_sf)
do elk come within 250 m of trails?
trails_and_roads_buff <- st_buffer(trails_and_roads2, 250)
library(mapview)
mapview(trails_and_roads2)+
mapview(trails_and_roads_buff, col.regions="red")+
mapview(subset(elk_combine, Used==TRUE))